Numerical Renormalization Group for Bosonic Systems 
and Application to the Subohmic Spin-Boson Model 
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We describe the generalization of Wilson's Numerical Renormalization Group method to quantum 
impurity models with a bosonic bath, providing a general non-perturbative approach to bosonic 
impurity models which can access exponentially small energies and temperatures. As an application, 
we consider the spin-boson model, describing a two-level system coupled to a bosonic bath with 
power- law spectral density, J(ui) oc uj s . We find clear evidence for a line of continuous quantum phase 
transitions for subohmic bath exponents < s < 1; the line terminates in the well-known Kosterlitz- 
Thouless transition at s — 1. Contact is made with results from perturbative renormalization group, 
and various other applications are outlined. 



The Numerical Renormalization Group method 
(NRG) developed by Wilson 

is a powerful tool for 
the investigation of the Kondo model and its generaliza- 
tions I n these models, a (possibly complex) 
impurity, such as a localized spin, couples to a fermionic 
bath. In the case of a spin-i impurity coupled antifer- 
romagnetically to a metallic bath, the impurity spin is 
screened below a characteristic scale Xk, the Kondo tem- 
perature |3J. The strength of the NRG lies in its non- 
perturbative nature and the ability to resolve arbitrarily 
small energies Q. A variety of thermodynamic and dy- 
namic quantities can be calculated for a large number of 
impurity models in the whole parameter space 0, 0] . 

There is, however, a very important class of models 
for which the NRG method has not yet been developed: 
models with a coupling of the impurity to a bosonic bath 
. The intensively studied spin-boson model 0, H| be- 
longs to this class; its Hamiltonian is given by 



H 



A 



<? z + knaiffli + \ (fflj + a\) ■ (i) 



Here the Pauli-matrices <7j describe a spin, i.e., a generic 
two-level system, which is linearly coupled to a bath of 
harmonic oscillators, with creation (annihilation) opera- 
tors a\ (a,i). The bare tunneling amplitude between the 
two spin states | f) and | |) is given by A, and e is an 
additional bias. The u>i are the oscillator frequencies and 
Ai the coupling strengths between the oscillators and the 
local spin. The coupling between spin and bosonic bath 
is completely specified by the bath spectral function 



J (ui) = ir X^5 (lo — LOi 



(2) 



Of particular interest are power-law spectra 



J[w) = 2iracjl s lo s , < w < lo c , s > -1 (3) 

where the dimensionless parameter a characterizes the 
dissipation strength, and lo c is a cutoff energy. The value 
s = 1 corresponds to the case of ohmic dissipation. 



The spin-boson model is a generic model describing 
quantum dissipation; it has been discussed in the con- 
text of a great variety of physical problems 0, @ rang- 
ing from the effect of friction on the electron transfer in 
biomolecules |9( to the description of the quantum entan- 
glement between a qubit and its environment |lft ITU Il2l | . 

Considering the wealth of applications, the question 
arises whether Wilson's NRG method can be exploited 
for this class of models; and it is the purpose of this paper 
to show that this is indeed the case. What we have in 
mind here is the direct mapping of models like Q to 
a semi-infinite chain form typical for the NRG |]|. As 
described below, bosonic operators constitute the sites 
of the chain, and the Hamiltonian is solved by iterative 
numerical diagonalization ^jjj . This approach is different 
from previous NRG calculations in R.efs. Hflll3. where the 
mapping of the spin-boson model Q to the anisotropic 
fermionic Kondo model was employed - such a mapping 
is restricted to the ohmic case J(u>) oc lo. 

Bosonic NRG. Let us now describe the generalization 
of the NRG method to a bosonic bath with a continu- 
ous spectrum. The strategy is similar to the one used 
for the Kondo or single- impurity Anderson model 0, • 
There are, however, important differences which we out- 
line here; a more detailed discussion will appear else- 
where. Here we present explicit equations for the spin- 
boson model p|l: the generalization to other impurity 
models or multiple bosonic baths is straightforward. 

We start from the following form of the model QJ : 



H = H loc + J deg(e)4a e + Y J deh(e){a e +4) (4) 



with H\ oc = —A<j x /2 + ea z /2. In this model, g(e) 
characterizes the dispersion of a bosonic bath in a one- 
dimensional representation, with upper cutoff 1 for e. 
The coupling between the spin and the bosonic bath is 
given by h(e). These two energy-dependent functions are 



2 



related to the spectral function J(co) via 

[e(x)\ (x€[0,w c 



b) 



-J{x) 

7T 



de(x) 2 



dx 



(5) 



where s(x) is the inverse function of g(x), g[e{x)\ = x. 
As discussed in Ref. for the Anderson model, Eq. © 
does not uniquely determine g(x) and h{x) 1 and a specific 
choice for h(x) is used to simplify the calculations. 

The NRG procedure starts by dividing the energy in- 
terval [0, 1] into intervals [A~(™ +1 ) , A~™] (n = 0, 1,2, . . .). 
An orthonormal set of functions ip n p(s) oc e l "™ pe is intro- 
duced for each interval so that the operators a £ can be 
represented in this basis. Choosing h{e) as constant in 
each interval and dropping the (p^O)-components as 
in 0, 0] , the Hamiltonian of the spin-boson model then 
takes the following form: 

oc oo 

H = #loc + ^ ?« a « a « + 53 7,1 ( a " + a ") ' 

■n— n V _,_ n 



n==0 

A _n a> e /-A~ 
2 



dxJ(a;)x, 7„ = / dxJ(x). (6) 



The transformation to a semi-infinite chain form yields 



'sem — A lloc 



lh, n — //|,„ + \J^~^ + ^o) 

n+l + 6 n+ i&n 

(7) 



n==0 



with rjo = J dxJ(x). The spin now couples to the first 
site of the bosonic chain only, and the remaining part of 
the chain is characterized by on-site energies e n and hop- 
ping parameters t n . The parameters t)q, e n and t n can 
be calculated numerically from a given spectral function 
J(w) 0|- Note that here the spectrum is restricted to 
positive frequencies; this results in hopping matrix ele- 
ments falling off as t n oc A~™ (which allows to work with 
A = 2 keeping a relatively small number of states), in 
contrast to the fermionic case where the discretization 
is performed for both negative and positive energies, in 
this case t n oc A~™/ 2 . The on-site energies also fall off 
as e n oc A~ n so that a fixed but s-dependent ratio t n /e n 
emerges for large n, where s is the bath exponent in |J3}. 

The Hamiltonian J7J is solved by iterative numerical 
diagonalization 0,0- At each step, one bosonic site of 
the chain is added. The infinite bosonic Hilbert space 
has to be cut off, by restricting the basis of each new 
bosonic site to a finite number of states Nj,. After diago- 
nalizing the enhanced cluster, the N s lowest lying many- 
particle states are kept, and the procedure is repeated 
|16| . The calculation of static and dynamic observables 
can be done in analogy to the fermionic NRG. In general, 
and as known from the fermionic case, the accuracy of 
the cutoff procedure has to be tested for each application, 
and we will show results below. 
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FIG. 1: a) Phase diagram for the transition between delocal- 
ized (a < a c ) and localized phases (a > a c ) of the spin-boson 
model Q for bias e = and various values of A, deduced from 
the NRG flow, b) A-dependence of the critical coupling a c 
for various values of the bath exponent s. The dashed lines 
are guides to the eye, the solid lines are fits to Eq. I1U1 using 
the A < 10~ 7 points only. For s close to 1 the asymptotic 
regime is only reached for very small A. NRG parameters 
here and in Figs. are A = 2, N b = 8, and iV s = 100. 



Application to the spin-boson model. To investigate 
the feasibility of the bosonic NRG, we performed exten- 
sive calculations for the spin-boson model with bath ex- 
ponents < s < 1, bias e = 0, and lo c = 1. In the 
ohmic case s = l, it is known that a Kosterlitz-Thouless 
quantum transition separates a localizedphase at a > a c 
from a delocalized phase at a < a c 0, ■ In the local- 
ized regime, the tunnel splitting between the two levels 
renormalizes to zero, whereas it is finite in the delocalized 
phase. For A -C lj c the transition occurs at a c = 1. 

The subohmic case 0,0] is less clear. For A/uj c — > 
the system is localized for any non-zero coupling, how- 
ever, the behavior at finite A was not discussed in 
Refs.0Jp For large A a delocalized phase was argued to 
exist ^ItTIi^I , and Ref. proposed a first-order transi- 
tion scenario. In the following, we shall resolve this issue 
and show that a continuous transition with associated 
critical behavior occurs for all < s < 1. 

Notably, the spin-boson model can be mapped onto a 
one-dimensional Ising model with long-range couplings 
falling off as r _s_1 ; the localized phase of the spin-boson 
model then corresponds to the ordered phase of the Ising 
magnet j?]. As shown by Dyson [Tg| . this Ising model 
features a transition for < s < 1, but the results for 
s < 1 have not been systematically carried over to the 
spin-boson model so far. 

Our NRG calculations provide clear evidence for a 
phase transition in the spin-boson model for all < s < 
1, which is continuous for < s < 1 and of Kosterlitz- 
Thouless type for s = 1. The numerical results are sum- 
marized in Fig. which shows the phase boundaries 
determined from the NRG flow for fixed NRG parame- 
ters A = 2, N b = 8, and N s = 100 [Tj|. (No transition 
occurs for s > 1: the system is always delocalized.) As 



FIG. 2: A dependence of the critical coupling a c at s — 0.9, 
A = 1CT 3 , for various NRG parameters N 3 and N b [l(J. The 
dashed line is a linear fit to the N 3 = 120, Nb — 8 data in the 
range 2 < A < 3. 



displayed in Fig. the critical coupling a c closely fol- 
lows a power law as function of the bare tunnel splitting, 
a c oc A x for small A, with an s-dependent exponent x. 
Our data are consistent with x = 1 — s, see below. 

In the ohmic case, s = 1, the critical a c approaches a fi- 
nite value as A — > 0. For the quoted NRG parameters we 
find a c s» 1.18, being slightly larger than the established 
value a c (s = 1, A — ► 0) = 1. This deviation is solely due 
to the NRG discretization; calculations with different A 
show that in the limit A —> 1 we recover a c = 1. The 
general behavior is illustrated in Fig. which shows a c 
for fixed A and s = 0.9. Keeping A fixed, we observe a 
rapid convergence of a c with increasing Nb and N s . As 
expected from the iterative diagonalization scheme, the 
values of Nb and N« n ecessary for convergence increase 
with decreasing A 16]. The converged data for a c (A) 
show a linear A dependence in the range 1.8 < A < 3, 
with a deviation of about 15% at A = 2 from the ex- 
trapolated A — > 1 value. The same holds for the ohmic 
case (data not shown) and the extrapolation results in 
a c (s = 1, A = 1CT 4 , A — ► 1) = 0.99 ± 0.02; our data are 
consistent with the RG result a c = 1 + 0(A/oj c ) 0. 

The NRG flow of the many-particle levels of the Hamil- 
tonian, displayed in Fig. |3| can be used to analyze the 
low-temperature behavior. For all values of < s < 1, 
we can identify two stable fixed points, corresponding to 
the localized and delocalized phases of the impurity spin, 
and a third NRG fixed point, which is infrared unstable 
and corresponds to a critical fixed point. In contrast, for 
s = 1 we find (in addition to the delocalized fixed point) 
a line of fixed points for a > a c , and no critical fixed 
point, as expected for a Kosterlitz-Thouless transition. 

The energy scale T* , describing the crossover from the 
critical to a stable fixed point, is shown in Fig. 0] For 
< s < 1, T* is found to vary in a power-law fashion 
with the distance from criticality, T* cx \ot — a c \ l/z , where 
we have introduced the correlation length and dynami- 
cal exponents v and z - note that l/(vz) is nothing but 
the scaling dimension of the leading relevant operator 
at the critical fixed point. Fig. HJ, nicely shows that vz 





FIG. 3: Sample NRG flow diagrams for the subohmic (upper 
panel, s — 0.6) and the ohmic case (lower panel, s = 1), close 
to the phase transition (except for the lower right panel). 



is independent of A for fixed s, further supporting the 
existence of a continuous quantum phase transition with 
universal behavior. In the ohmic case s = 1, T* varies ex- 
ponentially with the distance from the critical coupling, 
InT* cx l/(a c — a), as expected (Fig. Et>). 

In Fig. we show the s-dependence of the exponent 
vz. We find a divergence for both s — > and s — > 1, 
consistent with a Kosterlitz-Thouless transition at s = 1 , 
and the system being always localized at s = 0. The 
s — > 1 divergence is in good agreement with the pertur- 
bative result ©, see below. 

The NRG algorithm can be used to compute a vari- 
ety of static and dynamic observables. As an example, 
we show C(cl>), being the Fourier transform of the sym- 
metrized autocorrelation function C(t) — i([a z (t), o~ z ]+), 
in Fig. [SJj for s = 0.6 and parameters in the delocalized 
phase close to the transition. We observe a crossover 
from C(uS) oc lo s at small frequencies, characteristic of 
the delocalized phase 0, |2(| , to a quantum critical be- 
havior with a power-law divergence at higher frequencies 
- this gives rise to a characteristic peak at u> ~ T* . 
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FIG. 4: Crossover scale T* in the vicinity of the phase tran- 
sition for a < a c , defined here through the first excited NRG 
level A" E n i(T*) =0.3. a) Subohmic case s = 0.8, with power 
law fits, b) Ohmic case s = 1, with exponential fits. 
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FIG. 5: a) NRG results for the critical exponent vz, char- 
acterizing the vanishing of the crossover energy T* near the 
critical point, as function of the bath exponent s. Inset: data 
near s = 1 together with the function 1/^/2(1 — s) + 0.5, 
Eq. ©. b) Spin autocorrelation function C(uj) for s = 0.6, 
A = 10 -2 , and a = 0.19 < a c close to the transition. 

Comparison to perturbative results. The partition 
function of the spin-boson model can be approximately 
represented as that of a one-dimensional Ising model with 
couplings falling off as in this picture, defects in 

the Ising system correspond to spin flips of the original 
spin along the imaginary time axis. A RG analysis of 
this Ising model has been performed by Kosterlitz [2l| . 
Carrying over these results to the spin-boson model, we 
arrive at the RG equations (see also Ref. 11): 

0(a) = -a(A 2 + s - 1) , 

0(A) = A(l-a), (8) 

valid for small A, where A = A/lu c is the dimensionless 
tunneling strength. The RG flow is sketched in Fig. 1 of 
Ref. El For s — 1 , these equations are equivalent to the 
ones known from the anisotropic Kondo model, and de- 
scribe a Kosterlitz-Thouless transition, with a fixed line 
A = 0, a > 1. For s < 1 there is an unstable fixed point 
at a = 1, A 2 = 1 — s; clearly it is perturbatively acces- 
sible for small values of (1 — s) only. The critical fixed 
point is characterized by 

vz = 1/^2(1 - s) + 0(1) . (9) 

For small a, A, Eq. JHJ) yields for the phase boundary 

a c cx A 1 " 5 for A«w c (10) 

valid for all < s < 1. The results © and (0 are in 
good agreement with our numerical data in Figs.JS^, and 
respectively. 

Conclusions. We have presented a generalization 
of Wilson's NRG to quantum impurity problems with 
bosonic baths. Applying this novel technique to the 
subohmic spin-boson model, we have found a line of 
continuous boundary quantum phase transitions for all 
< s < 1, with exponents varying as function of s. This 



line terminates in a Kosterlitz-Thouless transition point 
at s = 1. Near s = 1, our numerical results are in agree- 
ment with perturbative calculations. The existence of a 
transition for s < 1 implies that weakly damped coher- 
ent dynamics is possible for qubits coupled to a subohmic 
bath, provided that the initial splitting A is large. 

In close analogy to the fcrmionic NRG, our method can 
be easily applied to the calculation of dynamical quan- 
tities. Furthermore, generalizations to impurities with 
multiple bosonic baths or both fcrmionic and bosonic 
baths are possible. This will allow the study of large 
classes of impurity models, e.g., so-called Bose-Kondo 
and Bose-Fermi-Kondo models. 
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